CUDA Python を使ってみよう

この資料はCUDA kernels in pythonを元に作成しています。

CUDA Python について

numba は、CUDA Pythonと呼ぶPythonコードの制限されたサブセットを提供しています。CUDA Python はCUDA実行モデルに従ってCUDAカーネルとデバイス関数にコンパイルすることで、CUDA GPUプログラミングをサポートします。 Numbaで記述されたカーネルは、NumPy配列に直接アクセスできるようです。 NumPy配列はCPUとGPUの間で自動的に転送されます。

Python でCUDAカーネルを実装

まず、配列内の各値の平方根を計算する最初のCUDAカーネルを実装してみましょう。

code: python

import numpy as np

a = np.arange(4096,dtype=np.float32)

https://gyazo.com/93da36631c3b4114ef60c434038ef6f7

numbaの@vectorizeデコレータを使用して、GPU上ですべての要素の平方根を並列計算することができます。

code: python

import math

from numba import vectorize

@vectorize('float32(float32)', target='cuda')

def gpu_sqrt(x):

return math.sqrt(x)

gpu_sqrt(a)

https://gyazo.com/bd15e38a53890b18bef85d636956be62

このコードをCUDAカーネルで実装してみましょう。

code: python

from numba import cuda

@cuda.jit

def gpu_sqrt_kernel(x, out):

idx = cuda.grid(1)

outidx = math.sqrt(xidx)

入力は4096の値をもつ配列aであるため、GPUで4096のスレッドを使用します。

入力配列と出力配列は1次元であるため、スレッドの1次元グリッドを使用します。 cuda.grid(1)を呼び出すと、グリッド全体の現在のスレッドのインデックス(一意の値)が返されます。 4096スレッドの場合、idx の範囲は0〜4095になります。

次に、コードで、各スレッドが入力配列の1つの要素を処理して、出力配列に1つの要素を生成することがわかります。この要素は、スレッドインデックスによってスレッドごとに決定されます。

CUDAカーネルができたので、入力配列をGPUデバイスにコピーし、同じ形状のデバイス上に出力配列を作成して、最後にカーネルを起動します。

code: python

# 入力データをGPUに転送

d_a = cuda.to_device(a)

# デバイス上に出力データを作成

d_out = cuda.device_array_like(d_a)

# それぞれ128スレッドを含む32ブロックを使用

blocks_per_grid = 32

threads_per_block = 128

gpu_sqrt_kernelblocks_per_grid, threads_per_block(d_a, d_out)

# すべてのスレッドが完了するのを待つ

cuda.synchronize()

# 出力配列をホストシステムに転送して出力

print(d_out.copy_to_host())

このコードでは、4096のスレッドが32ブロックのグリッドに配置され、各ブロックには128のスレッドがあります。

グリッドあたりのブロック数またはブロックあたりのスレッド数を小さくすると、スレッド数も小さくなります。この場合は、一部の要素は処理されず、出力配列の最後にある対応するスロットはデフォルト値のゼロ(0)に設定されたままになります。

一方、スレッド数を大きくすると、すべてが正常に動作しているように見えます。ただし、実際にはエラーが発生しています。エラーメッセージが常に表示されるとは限らないことに注意する必要があります。

最後に、cuda.synchronize()の呼び出しをコメントアウトすると、出力配列が部分的にしか満たされないか、まったく満たされないと予想するかもしれません。しかし、ホストへのコピーは暗黙的な同期を実行するため、実際には cuda.syncronize() を呼び出す必要はありません。

実行構成：ブロック数とスレッド数

この例では、4096スレッドを使用することにしました。これらのスレッドは、それぞれ128スレッドの32ブロックを含むグリッドとして配置されます。

次のような疑問があるかもしれません。

なぜグリッド内のブロックにスレッドを配置する必要があるか？

どうやって32と128というマジックナンバーを思いついたのか？

16と256（合計4096）や 64と64（合計4096）のような他の値を使用できるのか？

これらの疑問はパフォーマンスを気にする場合に重要なものです。

これらの疑問を明らかにするためには、GPUハードウェアについて理解する必要があります。

CUDA：GPUハードウェアについて学ぶ

cuda.detect()を使って使用しているGPUを調べてみましょう。

https://gyazo.com/77467d6d39e475ea03b47c5ffd32d3f7

この結果から、Tesla T4 GPU が使用されていることがわかります。

NVIDIAの資料によれば、Tesla T4 は Turing TU102GPUに基づいていることがわかります。

"GPUの概要とCUDAツールキット" で説明したように、T4ではデータ配列のサイズが約2万以上で最高のパフォーマンスを発揮できます。

はじめに例示したコードの配列aは4096でしたが、これは小さすぎるといえます。

ストライド

次の単純なカーネルは、入力配列から要素ひとつだけ処理します。

code: python

@cuda.jit

def gpu_atan(x, out):

idx = cuda.grid(1)

outidx = math.atan(xidx)

カーネルがデプロイされると、GPUは配列内の要素と同じ数のスレッドを作成しようとします。

配列が大きい場合は、多くのブロックが発生する可能性があります。

この場合、cuda.gridsize()で取得した数だけストライドさせ、ループを入力配列のいくつかの要素ことに処理するようにします。

code: python

@cuda.jit

def gpu_atan_stride(x, out):

start = cuda.grid(1)

stride = cuda.gridsize(1)

for i in range(start, x.shape0, stride):

outi = math.atan(xi)

このようにして、特定のスレッドが複数の要素を処理し、スレッド数が制御されます。スレッドは調整された方法で処理を続け、GPUはスレッドの作成とスケジューリングに時間を浪費しません。

次のような小さな例を考えてみましょう。

サイズ8の入力データ配列

それぞれ4スレッドのブロック

ストライドせずに、2つのブロックを使用する必要があるため、合計8つのスレッドがあり、それぞれが配列内の1つの要素を処理します。

code: python

a = np.arange(8, dtype=np.float32)

d_a = cuda.to_device(a)

d_out = cuda.device_array_like(d_a)

# 4スレッドの2ブロック

gpu_atan2,4(d_a, d_out)

print(d_out.copy_to_host())

ストライドを使用すると、4つのスレッドを持つ単一のブロックを使用できます。この場合、各スレッドは2つの要素を処理します。具体的には、上記のコードでは次のように処理されます。

startはスレッドインデックスであり、4つのスレッドでそれぞれ0、1、2、3。

strideはグリッドのサイズであり、グリッド内のスレッドの総数です。ここでは4。

x.shape[0]は、入力データ配列x のサイズであり、8。

code: python

# 4スレッドの1ブロック

gpu_atan_stride1,4(d_a, d_out)

print(d_out.copy_to_host())

各スレッドによって処理される要素は次のとおりです。

code: python

for thread_index in range(4):

print('thread', thread_index)

for j in range(thread_index, 8, 4):

print('\telement', j)

https://gyazo.com/38a8752afe05a60d6b8728352d17d183

グリッドが入力配列内のすべての要素を処理するために、次図のように移動していることを想像すると理解しやすいでしょう。

https://gyazo.com/16a5986295e1a0010b2d978eaa12728a

この図では、各色がスレッドに対応しています（たとえば、スレッド0は要素0と4を処理します）。

スレッドは非同期で進行することに注意してください。言い換えると、スレッド1が要素1を処理している間に、スレッド0が要素0と4で実行される可能性があります。

パフォーマンス分析

大きな配列の処理におけるストライドと実行構成パラメーターの影響を調べてみましょう。

ストライドの有無で２つのCUDAカーネルを実装します。

code: python

@cuda.jit

def gpu_atan(x, out):

idx = cuda.grid(1)

outidx = math.atan(xidx)

@cuda.jit

def gpu_atan_stride(x, out):

start = cuda.grid(1)

stride = cuda.gridsize(1)

for i in range(start, x.shape0, stride):

outi = math.atan(xi)

ここで、大きな配列を作成し、それをデバイスに出荷し、通常どおりデバイス上に出力配列を作成します。

code: python

import numpy as np

a = np.arange(256*1000000,dtype=np.float32)

# 入力データをGPUに転送

d_a = cuda.to_device(a)

# GPU上に出力配列を作成

d_out = cuda.device_array_like(d_a)

まず、GPUでこの配列を（シングルスレッドで）順次処理できる速度を見てみましょう。ここではストライドバージョンを使用します。これは、明らかに、単純バージョンでは1つのスレッドで1つの要素しか処理できないためです。

https://gyazo.com/ac9c547d304dd6a3324b87866c6b4937

GPU上の単一スレッドでこれらの2億5600万の値を処理するには、77秒かかりました。

並列処理がどのように役立つか見てみましょう。そのために、単純なルールに従って実行構成を選択し、ストライドしないCUDAカーネルを使用します。

https://gyazo.com/2308e1d68cf26a45b7bfc6d4f19dd566

並列バージョンは、予想どおりはるかに高速です。

それでは、ストライドバージョンによってパフォーマンスが向上するかどうかを試してみましょう。

https://gyazo.com/f89b69886106941956f1774922166637

ここでのパフォーマンスの改善率はそれほど重要ではありません。ばらつきがあるため、改善率を確認するためには、この2つのセルを数回再実行する必要がある場合があります。ストライドから得られるパフォーマンスの改善率は、入力配列のサイズと、カーネルによって実行される作業によって異なります。

多次元データセット：行列の乗算

これまで、スレッドの1Dグリッドを利用して、1次元配列を使用してきました。

次の1Dカーネルでは、cuda.grid(1)はグリッド内のスレッドの位置を識別する単一のインデックスを返し、cuda.gridsize(1)はスレッドのグリッドの長さを返します。

code: python

@cuda.jit

def gpu_atan_stride(x, out):

start = cuda.grid(1)

stride = cuda.gridsize(1)

for i in range(start, x.shape0, stride):

outi = math.atan(xi)

これらの関数は、2Dおよび3Dデータを処理するための簡単なインターフェイスも提供します。

スレッドの2Dグリッドをデプロイして2Dデータを処理する方法と、2Dデータをストライドする方法を説明することにしましょう。

シンプルな2Dカーネル

2Dデータを操作するための基本的なツールを説明するためみ。行列乗算を実装してみましょう。

code: python

@cuda.jit

def gpu_2d(out):

# 2Dでスレッド座標を取得

i1, i2 = cuda.grid(2)

outi1i2 = i1*10 + i2

見ての通り、このカーネルには入力はありません。その目的は、現在のスレッドの座標をエンコードして2D配列outに記録することだけです。

それでは、出力データ構造を作成しましょう。

https://gyazo.com/cd0de5124920aaa12aece5d02e84f98a

このデータ構造をGPUにコピーし、カーネルをデプロイして、結果を取り出します。

https://gyazo.com/a7e21f5258fb7e59fdab14ee4e10eccd

最初の2つの列はスレッドの最初のブロックによって処理され、最後の2つは2番目のブロックによって処理されました。

カーネルコードとこの出力を見ると、最初の次元（行）がi1にインデックス付けられ、2番目の次元（列）がi2にインデックス付けられていることがわかります。

1Dの場合と同様に、グリッドが出力データ構造全体をカバーしていることを確認する必要があります。たとえば、ブロック構造を考慮しない場合、最後の2つの列は未処理のままであり、実際には2D配列の「下」にある許可されていないメモリにアクセスします。

https://gyazo.com/cc4266fa2f65c20cc6cf352fcfc5c9e0

2Dストライド

2Dでのストライドは、1Dよりもそれほど難しくありません。繰り返しになりますが、非常に簡単な例を見てみましょう。まず、ストライドの実装が機能していることを確認します。

code: python

@cuda.jit

def gpu_2d_stride(out):

# 2Dでスレッド座標を取得

s1, s2 = cuda.grid(2)

# 2Dでグリッドサイズを取得

d1, d2 = cuda.gridsize(2)

for i1 in range(s1, out.shape0, d1):

for i2 in range(s2, out.shape1, d2):

outi1i2 = i1*10 + i2

デバイスの出力データをリセットし、(3,2)ブロックをひとつ展開します。したがって、グリッドは出力配列の半分のみをカバーします。ただし、ストライドは出力配列内のグリッドを移動するため、配列全体をカバーします。これを確認してみましょう。

https://gyazo.com/409afcb4926e5ac6e87fb45737e53595

ストライドしない場合と同じ結果が得られます。次に、各スレッドの開始点の指定を(i1, i2) から(s1, s2)に変更したカーネルでストライドの挙動を確認してみます。

code: python

@cuda.jit

def gpu_2d_start(out):

# 2Dでスレッド座標を取得

s1, s2 = cuda.grid(2)

# 2Dでグリッドサイズを取得

d1, d2 = cuda.gridsize(2)

for i1 in range(s1, out.shape0, d1):

for i2 in range(s2, out.shape1, d2):

# 以前のバージョンのカーネルとの違いに注意

outi1i2 = s1*10 + s2

https://gyazo.com/4f188ec92f549ceaa419657c6d8b42fa

同じ値の要素は同じスレッドで処理され、右側に1つストライドしていることがわかります。

2D配列を拡張すると、ストライドが両方向に発生することがわかります。

https://gyazo.com/99bd6feec83452d4d28583994d627e43

ここでは、(2,2)構造に配置された4つのストライドがあります。下部のストライドは不完全ですが、ループが出力配列の次元に制限されているため、カーネルはオーバーフローから自然に保護されます。

2D行列の乗算：理論

2D行列乗算カーネルを実装することは、2Dでのストライドを理解しているかを確認するための良い方法になります。

まず、理論について思い出してみましょう。サイズ(m, n)の2つの行列Aとサイズ(n, p)の行列Bは、行列Aの列の数が行列Bの行の数に等しいため、乗算することができます。積行列Cのサイズ(m, p)になります。そして、行列Cの係数の1つの値は、次の式で与えられます。

$ c_{ij} = \sum_{k=0}^{n} a_{ik} b_{kj}

線形代数を初めてみる場合、あるいは長い間使用していなかったあ場合、これは少し難しくに見えるかもしれません。それでは、コードの簡単な例を見てみましょう。 2つの行列を作成します。

https://gyazo.com/2ba2fca2a49b73d2078c5aa0f2aef1d6

これ、numpyを使用してCPUで乗算を行います。

https://gyazo.com/732987495feeebed63f34de59d000a7e

行列aのサイズは(2, 3)であり、行列bのサイズは(3, 4)なので、積行列cのサイズは(2, 4)になります。

たとえば、積行列cの1行目と2列目の係数である23について考えてみましょう。これはc01です。上記の和分方程式は、この係数は、すべてを合計する前に、行列aの最初の行の係数を取得し、それらに行列bの2番目の列の係数をそれぞれ乗算することによって得られることを示しています。

次のようになります：c01 = a00*b01 + a01*b11 + a02*b21 = 0*1 + 1*5 + 2*9 = 23

2D行列乗算カーネル

最初のステップとして、ストライドせずに行列の乗算を行うカーネルを実装します。つまり、出力マトリックスの要素ごとに1つのスレッドが必要になります。カーネルは次のとおりです。

code: python

@cuda.jit

def multiply(a, b, c):

i1, i2 = cuda.grid(2)

the_sum = 0

for k in range(b.shape0): # または a.shape1 同じ

the_sum += ai1k*bki2

ci1, i2 = the_sum

行列aとbをGPUに転送し、出力行列cを作成します。

https://gyazo.com/36fbc762992110b418177e3fda3492ba

次に、カーネルをデプロイします。出力マトリックスの形状に一致する、2行4列の2D配列として配置された、8つのスレッドを持つ単一のブロックを使用します。

https://gyazo.com/c55fb18ead79d1b00ecfb70eedf5be10

いいですね、numpyと同じ結果が得られています。

ストライドの実装は実際には大したことではありません。

少し注意が必要なのことは、(a.shape[0], b.shape[1])の形状を持つ出力行列をまたぐ必要があることを覚えておくことだけです。ストライドの場合、入力は関係ありません。

code: python

@cuda.jit

def multiply_stride(a, b, c):

s1, s2 = cuda.grid(2)

d1, d2 = cuda.gridsize(2)

for i1 in range(s1, a.shape0, d1):

for i2 in range(s2, b.shape1, d2):

the_sum = 0

for k in range(b.shape0): # a.shape1 と同じ

the_sum += ai1k*bki2

ci1, i2 = the_sum

https://gyazo.com/34c28163643161d05c67852c129d269d

その他のCUDAライブラリ

CUDA Python では次のCUDAライブラリも利用することができます。

cuDNN：ディープニューラルネットワーク用のプリミティブのGPUアクセラレーションライブラリ。バックプロパゲーション層や畳み込み層を自分でコーディングしなくてもよくなる。

cuBLAS：GPU対応の線形代数ライブラリ

cuSPARSE：GPU対応のスパース行列（疎行列)のための線形代数演算

cuRAND：GPU対応の乱数生成

cuFFT：GPU対応の高速フーリエ変換

CUDA Pythonの詳細

実行モデル

CUDA Pythonは、CUDAのSIMT(Single Instruction Multiple Thread)モデルに直接マップします。各命令は、複数のスレッドによって暗黙的に並列に実行されます。この実行モデルでは、複数のスレッドが同じタスクを実行する必要がないため、配列式はあまり役に立ちません。代わりに、スレッドが協調してタスクを実行するようにします。

ホスト説明: SIMT（Single Instruction Multiple Thread）

NVIDIAは2006年11月に発表したG80 GPUから、一般的なプロセッサとは異なる SIMTアーキテクチャを採用しています。SIMTは、汎用プロセッサにおけるSIMD(single instruction, multiple data)とハードウェアマルチスレッドの中間的な構造をしており、1つの命令デコーダで複数の実行ユニットを制御しますが、SIMDとは違い各処理単位が独立したarchitectural stateを持ちます。

CUDA PythonでサポートされているPython機能について説明しましょう。

構文

次のPythonの構文はサポートされていません。

例外処理 (try .. except, try .. finally)

コンテキストマネージャー(with文)

内包表記 (リスト、辞書、セット、ジェネレータの内包表記)

ジェネレーター（任意のyieldステートメント）

次のPythonの構文はサポートされています。

raise文

assert文

debug=Trueが numba.cuda.jit()デコレータに渡された場合にのみ効果があります。

CUDA C / C ++のassertキーワードの動作に似ています。これは、デバイスのデバッグをオンにしてコンパイルしない限り無視されます。

文字列、整数、および浮動小数点数のprint()はサポートされています

print()は非同期操作です。

カーネルの起動後にすべての出力が確実に印刷されるようにするには、numba.cuda.synchronize() を呼び出して同期する必要があります。

同期処理を省略できますが、カーネルからの出力のタイミングは不定になります。

その後に実行されるドライバー操作中に表示される場合がある

プログラムの実行が完了する前に表示されない場合がある

組み込み型

次の組み込み型のサポートは、CPU nopythonモードから継承されます。

int

float

complex

bool

tuple

None

組み込み関数

次の組み込み関数がサポートされています。

abs()

bool()

complex()

enumerate()

float()

int()：引数は１つだけの呼び出しのみ

len()

min(): 引数が複数の呼び出しのみ

max(): 引数が複数の呼び出しのみ

range()

round()

zip()

標準ライブラリモジュール

cmathモジュール

cmathモジュールの次の機能がサポートされています。

cmath.acos()

cmath.acosh()

cmath.asin()

cmath.asinh()

cmath.atan()

cmath.atanh()

cmath.cos()

cmath.cosh()

cmath.exp()

cmath.isfinite()

cmath.isinf()

cmath.isnan()

cmath.log()

cmath.log10()

cmath.phase()

cmath.polar()

cmath.rect()

cmath.sin()

cmath.sinh()

cmath.sqrt()

cmath.tan()

cmath.tanh()

mathモジュール]

math モジュールの次の関数がサポートされています。

math.acos()

math.asin()

math.atan()

math.acosh()

math.asinh()

math.atanh()

math.cos()

math.sin()

math.tan()

math.hypot()

math.cosh()

math.sinh()

math.tanh()

math.atan2()

math.erf()

math.erfc()

math.exp()

math.expm1()

math.fabs()

math.gamma()

math.lgamma()

math.log()

math.log10()

math.log1p()

math.sqrt()

math.pow()

math.ceil()

math.floor()

math.copysign()

math.fmod()

math.modf()

math.isnan()

math.isinf()

operatorモジュール

operatorモジュールの次の機能がサポートされています。

operator.add()

operator.and_()

operator.eq()

operator.floordiv()

operator.ge()

operator.gt()

operator.iadd()

operator.iand()

operator.ifloordiv()

operator.ilshift()

operator.imod()

operator.imul()

operator.invert()

operator.ior()

operator.ipow()

operator.irshift()

operator.isub()

operator.itruediv()

operator.ixor()

operator.le()

operator.lshift()

operator.lt()

operator.mod()

operator.mul()

operator.ne()

operator.neg()

operator.not_()

operator.or_()

operator.pos()

operator.pow()

operator.rshift()

operator.sub()

operator.truediv()

operator.xor()

Numpyサポート

CUDAプログラミングモデルのため、カーネル内の動的メモリ割り当ては非効率的であり、多くの場合必要ありません。 Numbaはメモリ割り当て機能を許可していません。これにより、多数のNumPyAPIが無効になります。最高のパフォーマンスを得るには、各スレッドが一度に1つの要素を処理するようにコードを記述する必要があります。

サポートされているNumPyの機能

ndarrayアトリビュートへのアクセス: .shape、.strides、.ndim、.sizeなど

mathモジュールに同等のものがあるスカラーufunc：np.sin(x[0])、xは1D配列

インデックス作成とスライス操作

サポートされていないNumPyの機能

ndarray作成API

ndarryのメソッド

新しいndarrayを返す関数

アトミック操作

Numbaは、CUDAでサポートされているいくつかのアトミック操作へのアクセスを提供します。現在実装されているものは次のとおりです。

class numba.cuda.atomic

アトミック操作の名前空間

class add(ary, idx, val)

アトミックary[idx] + = val を実行します。 int32、float32、およびfloat64オペランドでのみサポートされます。

アトミックにロードされているかのように、インデックス位置の古い値を返します。

class compare_and_swap(ary, old, val)

現在の値がoldと一致する場合は、条件付きでvalを１次元配列の最初の要素に割り当てます。

アトミックにロードされているかのように現在の値を返します。

class max(ary, idx ,val)

アトミックary[idx] = max(ary[idx], val) を実行します。

int32、int64、uint32、uint64、float32、float64オペランドでのみサポートされます。

アトミックにロードされているかのように、インデックスの場所にある古い値を返します。

class min(ary, idx, val)

アトミックary[idx] = min(ary[idx], val)を実行します。

int32、int64、uint32、uint64、float32、float64オペランドでのみサポートされます。

アトミックにロードされているかのように、インデックスの場所にある古い値を返します。

次のコードは、numba.cuda.atomic.maxを使用して配列の最大値を見つける方法を示しています。この場合、これは最大値を見つける最も効率的な方法ではありませんが、例として役立つことに注意してください。

code: pytohn

from numba import cuda

import numpy as np

@cuda.jit

def max_example(result, values):

"""Find the maximum value in values and store in result0"""

tid = cuda.threadIdx.x

bid = cuda.blockIdx.x

bdim = cuda.blockDim.x

i = (bid * bdim) + tid

cuda.atomic.max(result, 0, valuesi)

arr = np.random.rand(16384)

result = np.zeros(1, dtype=np.float64)

max_example256,64(result, arr)

print(result0) # Found using cuda.atomic.max

print(max(arr)) # Print max(arr) for comparison (should be equal!)

複数次元配列は、インデックスにintのタプルを使用することでサポートされます。

code: pyton

@cuda.jit

def max_example_3d(result, values):

"""

Find the maximum value in values and store in result0.

Both result and values are 3d arrays.

"""

i, j, k = cuda.grid(3)

# Atomically store to result0,1,2 from valuesi, j, k

cuda.atomic.max(result, (0, 1, 2), valuesi, j, k)

arr = np.random.rand(1000).reshape(10,10,10)

result = np.zeros((3, 3, 3), dtype=np.float64)

max_example_3d(2, 2, 2), (5, 5, 5)(result, arr)

print(result0, 1, 2, '==', np.max(arr))

乱数の生成

Numbaは、GPUで実行できる乱数生成アルゴリズムを提供します。ただし、NVIDIAがcuRANDを実装する方法に関する技術的な問題から、NumbaのGPU乱数ジェネレーターはcuRANDに基づいていません。代わりに、NumbaのGPURNGはxoroshiro128+アルゴリズムを実装したものです。 xoroshiro128+アルゴリズムの周期は$ 2^{128}-1 であり、これはcuRANDでデフォルトで使用されるXORWOWアルゴリズムの周期よりも短いですが、xoroshiro128+は乱数ジェネレーター品質のBigCrushテストに合格しています。

GPURNGを使用する場合は、各スレッドに独自のRNG状態があり、重複しないシーケンスを生成するように初期化されていることを確認することが重要です。 numba.cuda.randomモジュールは、これを行うためのホスト関数と、均一または正規分布の乱数を取得するためのCUDAデバイス関数を提供します。

numba.cuda.random.create_xoroshiro128p_states(n, seed, subsequence_start=0, stream=0)

n個の乱数ジェネレーター用に初期化された新しいデバイス配列を返します。

これにより、RNG状態が初期化され、配列内の各状態が、メインシーケンス内で互いに$ 2^{64}ステップ離れたサブシーケンスに対応するようになります。したがって、CUDAスレッドが$ 2^{64}を超える乱数を要求しない限り、この関数によって生成されるすべてのRNG状態は独立していることが保証されます。

n: int 作成するRNG状態の数

seed: uint64 ジェネレーターのリストの開始シード

subsequence_start: uint64 最初のRNG状態を$ 2^{64}ステップの倍数だけ進める

stream: CUDA Steam　初期化カーネルを実行するストリーム

numba.cuda.random.init_xoroshiro128p_states(states, seed, subsequence_start=0, stream=0)

並列ジェネレーターのGPUでRNG状態を初期化します。

subsequence_startパラメーターを使用して、最初のRNG状態を$ 2^{64}ステップの倍数だけ進めることができます。

state: 1D DeviceNDArray, dtype = xoroshiro128p_dtype　RNG状態の配列

seed: uint64 ジェネレーターのリストの開始シード

numba.cuda.random.xoroshiro128p_uniform_float32(states, index)

range(0.0, 1.0) のfloat32を返し、states[index] を進めます。

state: 1D array, dtype = xoroshiro128p_dtype RNG状態の配列

index: int64] 更新する状態のオフセット

戻り値の型: float32

numba.cuda.random.xoroshiro128p_uniform_float64(states, index)

range(0.0 ,1.0) のfloat64を返し、states[index] を進めます。

state: 1D array, dtype = xoroshiro128p_dtype RNG状態の配列

index: int64 更新する状態のオフセット

戻り値の型: float64

numba.cuda.random.xoroshiro128p_normal_float32(states, index)

正規分布のfloat32を返し、states[index] を進めます。

戻り値は、Box-Muller変換を使用し、mean=0 およびsigma=1 のガウス関数から取得されます。

これにより、RNGシーケンスが2ステップ進みます。

state: 1D array, dtype = xoroshiro128p_dtype RNG状態の配列

index: int64 更新する状態のオフセット

戻り値の型 float32

numba.cuda.random.xoroshiro128p_normal_float64(states, index)

正規分布のfloat32を返し、states[index] を進めます。

戻り値は、Box-Muller変換を使用して、mean=0 および sigma=1 のガウス関数から取得されます。これにより、RNGシーケンスが2ステップ進みます。

states: 1D array, dtype = xoroshiro128p_dtype RNG状態の配列

index: int64 更新する状態のオフセット

戻り値の型 float64

code: python

from __future__ import print_function, absolute_import

from numba import cuda

from numba.cuda.random import create_xoroshiro128p_states, xoroshiro128p_uniform_float32

import numpy as np

@cuda.jit

def compute_pi(rng_states, iterations, out):

"""値の最大値を見つけresult0に保存"""

thread_id = cuda.grid(1)

# ランダムな点(x ,y)を描画し、単位円内にある点の比率から円周率を計算

inside = 0

for i in range(iterations):

x = xoroshiro128p_uniform_float32(rng_states, thread_id)

y = xoroshiro128p_uniform_float32(rng_states, thread_id)

if x**2 + y**2 <= 1.0:

inside += 1

outthread_id = 4.0 * inside / iterations

threads_per_block = 64

blocks = 24

rng_states = create_xoroshiro128p_states(threads_per_block * blocks, seed=1)

out = np.zeros(threads_per_block * blocks, dtype=np.float32)

compute_piblocks, threads_per_block(rng_states, 10000, out)

print('pi:', out.mean())

デバイス管理

マルチGPUマシンの場合、ユーザーは使用するGPUを選択することができます。デフォルトでは、CUDAドライバーはデバイス0として最速のGPUを選択します。これは、Numbaが使用するデフォルトのデバイスです。

ここでの説明は、CUDA対応の複数のGPUが搭載されたシステムでのみ必要なことです。

デバイスの選択

必要な場合は、CUDA機能を使用する前にデバイスの選択を行う必要があります。

code: python

from numba import cuda

cuda.select_device(0)

デバイスは次の方法でクローズすることができます。

code: python

cuda.close()

その後、ユーザーは別のデバイスで新しいコンテキストを作成できます。

code: python

cuda.select_device(1) # GPUが2つあるとして

numba.cuda.select_device(device_id)

選択したdevice_idの新しいCUDAコンテキストを作成します。 device_idは、デバイスの番号である必要があります（0から始まり、デバイスの順序はCUDAライブラリによって決定されます）。コンテキストは現在のスレッドに関連付けられています。 Numbaは現在、スレッドごとに1つのコンテキストのみを許可しています。

成功した場合、この関数はデバイスインスタンスを返します。

numba.cuda.close()

現在のスレッドのすべてのコンテキストを明示的に閉じます。

デバイスリスト

デバイスリストは、システム内のすべてのGPUのリストであり、インデックスを作成して、選択したGPUでの実行を保証するコンテキストマネージャーを取得することができます。

numba.cuda.gpus

numba.cuda.cudadrv.devices.gpus

numba.cuda.gpus は、_DeviceListクラスのインスタンスであり、そこから現在のGPUコンテキストを取得することもできます。

class numba.cuda.cudadrv.devices._DeviceList

プロパティ: current

アクティブなデバイスを返します。アクティブなデバイスがない場合はNoneを返します

行列の乗算

次のコードは、CUDAカーネルを使用した行列乗算の単純な実装例です。

code: python

@cuda.jit

def matmul(A, B, C):

""" C = A * B　の二乗行列乗算を実行

"""

i, j = cuda.grid(2)

if i < C.shape0 and j < C.shape1:

tmp = 0.

for k in range(A.shape1):

tmp += Ai, k * Bk, j

Ci, j = tmp

この実装は簡単で直感的ですが、同じマトリックス要素がデバイスメモリから複数回ロードされるため、パフォーマンスが低下します。一部のデバイスには透過的なデータキャッシュがある場合がありますが、入力全体を一度に保持するのに十分な大きさではない場合があるため、これは遅くなります。

ブロックアルゴリズムを使用してデバイスメモリへのアクセスを減らすと、もっと高速になります。 CUDAは、ブロック内のスレッドに高速共有メモリを提供して、タスクを協調的に計算します。

次のコードは、共有メモリを使用した正方行列乗算の高速バージョンを実装したものです。

code: python

from numba import cuda, float32

# ブロックごとのスレッドと共有メモリの使用量を制御する

# TPBxTPB 要素のブロックで計算する

TPB = 16

@cuda.jit

def fast_matmul(A, B, C):

# 共有メモリに配列を定義

# 配列のサイズと型は、コンパイル時に既知であることが必要

sA = cuda.shared.array(shape=(TPB, TPB), dtype=float32)

sB = cuda.shared.array(shape=(TPB, TPB), dtype=float32)

x, y = cuda.grid(2)

tx = cuda.threadIdx.x

ty = cuda.threadIdx.y

bpg = cuda.gridDim.x # グリッドあたりのブロック数

if x >= C.shape0 and y >= C.shape1:

# (x, y) が有効なC境界の外側にある場合は終了

return

# 各スレッドは、結果マトリックスの1つの要素を計算

# 内積はTPB長ベクトルの内積にチャンク化される

tmp = 0.

for i in range(bpg):

# 共有メモリにデータをプリロード

sAtx, ty = Ax, ty + i * TPB

sBtx, ty = Btx + i * TPB, y

# すべてのスレッドがプリロードを完了するまで待つ

cuda.syncthreads()

# 共有メモリ上の部分積を計算

for j in range(TPB):

tmp += sAtx, j * sBj, ty

# すべてのスレッドが計算を終了するまで待つ

cuda.syncthreads()

Cx, y = tmp

共有メモリは限られたリソースであるため、コードは入力配列から一度に小さなブロックをプリロードします。次に、syncthreads() を呼び出して、すべてのスレッドがプリロードを完了するまで待機し、共有メモリで計算を実行します。計算後に再度同期して、次のループ反復でデータを上書きする前に、すべてのスレッドが共有メモリ内のデータで終了したことを確認します。

CUDAシミュレーターを使用したCUDAPythonのデバッグ

Numbaには、Pythonインタープリターといくつかの追加のPythonコードを使用して、CUDA Pythonのほとんどのセマンティクスを実装するCUDAシミュレーターが含まれています。これは、コードにprint()を追加するか、デバッガーを使用して個々のスレッドの実行をステップ実行することにより、CUDA Pythonコードをデバッグするために使用できます。

カーネルの実行は、シミュレーターによって一度に1ブロックずつ実行されます。ブロック内のスレッドごとに1つのスレッドが生成され、これらのスレッドの実行のスケジューリングはオペレーティングシステムに任されています。

シミュレーターの使用

シミュレーターは、環境変数NUMBA_ENABLE_CUDASIMを1に設定することで有効になります。その後、CUDA Pythonコードを通常どおり実行できます。カーネル内でデバッガーを使用する最も簡単な方法は、単一のスレッドのみを停止することです。そうしないと、デバッガーとの対話を処理するのが困難になります。たとえば、以下のカーネルはスレッド<<<(3,0,0), (1, 0, 0)>>>で停止します。

code: python

@cuda.jit

def vec_add(A, B, out):

x = cuda.threadIdx.x

bx = cuda.blockIdx.x

bdx = cuda.blockDim.x

if x == 1 and bx == 3:

from pdb import set_trace; set_trace()

i = bx * bdx + x

outi = Ai + Bi

1次元グリッドと1次元ブロックで呼び出された場合。

サポートされている機能

シミュレーターは、実際のGPUでの実行を可能な限り完全にシミュレーションすることを目的としています。特に、以下がサポートされています。

不可分操作

コンスタントメモリ

ローカルメモリ

共有メモリ

シミュレータはソース行情報を使用してスレッド間の共有メモリの割り当てを追跡するため、共有メモリ配列の宣言は別々のソース行に行う必要があります。

syncthreads()がサポートされていますが、異なるスレッドが異なるsyncthreads()呼び出しを入力した場合、起動は失敗しませんが、予期しない動作が発生します。シミュレータの将来のバージョンでは、この状態が検出される可能性があります。

ストリームAPIはサポートされていますが、実際のデバイスとは異なり、すべての操作は順次同期的に行われます。したがって、ストリームでの同期は操作できません。

イベントAPIもサポートされていますが、意味のあるタイミング情報は提供されません。

GPUとの間のデータ転送-特に、device_array()とdevice_array_like()を使用して配列オブジェクトを作成します。ピン留めされたメモリpinned()およびpinned_array()のAPIもサポートされていますが、ピン留めは行われません。

GPUコンテキストのリスト（cuda.gpusおよびcuda.cudadrv.devices.gpus）のドライバーAPI実装がサポートされており、単一のGPUコンテキストを報告します。このコンテキストは、実際のコンテキストと同じようにクローズしてリセットできます。

detect()関数がサポートされており、SIMULATORと呼ばれる1つのデバイスを報告します。

シミュレータの制限

シミュレータのいくつかの制限があります。

型チェック/型推論は行わない

JIT関数の引数の型が正しくない場合やまたはローカル変数の型の指定が正しくない場合、シミュレーターによって検出されません。

1つのGPUのみがシミュレートされます。

単一のGPUへのマルチスレッドアクセスはサポートされていません。

ほとんどのドライバーAPIは実装されていません。

PTXコードをCUDAPython関数とリンクすることはできません。

ワープおよびワープレベルの操作はまだ実装されていません。

明らかに実際のデバイスの速度よりもはるかに低速になる

シミュレータでのデバッグを扱いやすくするために、入力データのサイズとCUDAグリッドのサイズを小さくする必要がある場合があります。

GPUの縮約

CUDA GPUで縮約アルゴリズム(Reduction Algorithm)を作成するのは難しい場合があります。 Numbaは、単純な二項演算をリダクションカーネルに変換するための@reduceデコレータを提供します。

code: python

import numpy

from numba import cuda

@cuda.reduce

def sum_reduce(a, b):

return a + b

A = (numpy.arange(1234, dtype=numpy.float64)) + 1

expect = A.sum() # numpy sum reduction

got = sum_reduce(A) # cuda sum reduction

assert expect == got

ここで、ラムダ式を使用することができます：

code: python

sum_reduce = cuda.reduce(lambda a, b: a + b)

補足説明：縮約(REDUCTION)

並列計算においては，複数のスレッドやプロセスが保持する数値の総和を計算することを言い、

様々な並列計算で頻繁に用いられる、この Reductionはスレッド数やプロセス数に応じて、

処理時間が増えるため，高速化は重要な課題になります。

class Reduce

@reduceデコレータは、Reduceクラスのインスタンスを作成します。

class numba.cuda.Reduce(functor)

__call __(self, arr, size=None, res=None ,init=0, stream=0)

完全な削減を実行します。

arr ：ホストまたはデバイスアレイ。デバイス配列が指定されている場合、リダクションはインプレースで実行され、配列の値が上書きされます。ホストアレイが指定されている場合は、デバイスに自動的にコピーされます。

size：削減するarrの要素数を指定するオプションの整数。このパラメーターが指定されていない場合、配列全体が縮小されます。

res：削減結果を書き込むオプションのデバイス配列。結果は、この配列の最初の要素に書き込まれます。このパラメーターが指定されている場合、デバイスからホストへのリダクション出力の通信は行われません。

init：リダクションのオプションの初期値。タイプはarr.dtypeと一致する必要があります。

stream：削減を実行するオプションのCUDAストリーム。ストリームが指定されていない場合、デフォルトのストリーム0が使用されます。

戻り値

resが指定されている場合、Noneが返されます。それ以外の場合は、削減の結果が返されます。

__init __(self, functor)

指定されたバイナリ関数を使用して値を削減する削減オブジェクトを作成します。バイナリ関数は一度コンパイルされ、このオブジェクト内にキャッシュされます。このオブジェクトを存続させると、再コンパイルが防止されます。

binop：CUDAデバイスでの削減のための二項演算として使用されるCUDAデバイス関数としてコンパイルされる関数。内部的には、cuda.jit(device=True) を使用してコンパイルされます。

CUDA Ufuncs と Ufuncs

このページでは、CUDA Ufuncのようなオブジェクトについて説明します。

CUDAプログラムのプログラミングパターンをサポートするために、CUDAVectorizeとGUVectorizeは従来のufuncを生成できません。代わりに、ufuncのようなオブジェクトが返されます。このオブジェクトは類似したものですが、通常のNumPyufuncとは完全には互換性がありません。 CUDA ufuncは、デバイス内アレイ（すでにGPUデバイス上にある）を渡すためのサポートを追加して、PCI expressバス上のトラフィックを削減します。また、非同期モードで起動するためのstreamキーワード引数を受け入れます。

code: python

import math

from numba import vectorize, cuda

import numpy as np

@vectorize(['float32(float32, float32, float32)',

'float64(float64, float64, float64)'],

target='cuda')

def cu_discriminant(a, b, c):

return math.sqrt(b ** 2 - 4 * a * c)

N = 10000

dtype = np.float32

# prepare the input

A = np.array(np.random.sample(N), dtype=dtype)

B = np.array(np.random.sample(N) + 10, dtype=dtype)

C = np.array(np.random.sample(N), dtype=dtype)

D = cu_discriminant(A, B, C)

print(D)

デバイス関数の呼び出し

すべてのCUDAufuncカーネルには、他のCUDAデバイス関数を呼び出す機能があります。

code: python

from numba import vectorize, cuda

# デバイス機能を定義する

@cuda.jit('float32(float32, float32, float32)', device=True, inline=True)

def cu_device_fn(x, y, z):

return x ** y / z

# デバイス関数を呼び出すufuncを定義します

@vectorize('float32(float32, float32, float32)', target='cuda')

def cu_ufunc(x, y, z):

return cu_device_fn(x, y, z)

一般化されたCUDA ufuncs

一般化されたufuncは、CUDA ufunc機能と同様に、CUDAを使用してGPUで実行できます。これは、次のように実行できます。

code: python

from numba import guvectorize

@guvectorize(['void(float32:,:, float32:,:, float32:,:)'],

'(m,n),(n,p)->(m,p)', target='cuda')

def matmulcore(A, B, C):

# ...

gufuncカーネルがGPUのリソースを使いすぎて、カーネルの起動が失敗する場合があります。ユーザーは、コンパイルされたgufuncオブジェクトにmax_blocksizeアトリビュートを設定することにより、スレッドブロックの最大サイズを明示的に制御できます。

code: python

from numba import guvectorize

@guvectorize(..., target='cuda')

def very_complex_kernel(A, B, C):

# ...

very_complex_kernel.max_blocksize = 32 # ブロックあたり32スレッドに制限

CUDAメモリの共有

プロセス間の共有

注意：この機能はLinuxのみに制限されています。

デバイスアレイを別のプロセスにエクスポートする

デバイスアレイは、CUDA IPC APIを使用して、同じマシン内の別のプロセスと共有できます。これを行うには、デバイス配列でget_ipc_handle()メソッドを使用して、別のプロセスに転送できるIpcArrayHandleオブジェクトを取得します。

DeviceNDArray.get_ipc_handle(self)

シリアル化して別のプロセスに転送して、ローカル割り当てを共有しても安全なIpcArrayHandleオブジェクトを返します。

class numba.cuda.cudadrv.devicearray.IpcArrayHandle(ipc_handle, array_desc)

GPU割り当てを共有するために、シリアル化して同じマシン内の別のプロセスに転送できるIPC配列ハンドル。

転送先のプロセスで、open()メソッドを使用して、転送元のプロセスからの割り当てを共有する新しいDeviceNDArrayオブジェクトを作成します。リソースを解放するには、close()メソッドを呼び出します。その後、宛先は共有配列オブジェクトを使用できなくなります。

このオブジェクトは、open() と close() メソッドを自動的に呼び出すコンテキストマネージャーインターフェイスを実装します。

code: python

with the_ipc_array_handle as ipc_array:

# ここでipc_arrayを通常のGPU配列オブジェクトとして使用する

some_code(ipc_array)

# ipc_array はこの時点で廃棄される

close(self)

GPU arrayへのIPCハンドルを閉じます。

open(self)

転送元のプロセスからの割り当てを共有する新しいDeviceNDArrayを返します。転送元のプロセスでは使用しないでください。

別のプロセスからIPCメモリをインポート

次の関数は、デバイス配列として別のプロセスからIPCハンドルを開くために使用されます。

cuda.open_ipc_arra(handle, shape, dtype, strides=None, offset=0)

バイトシーケンス（バイト、intのタプルなど）として表されるIPCハンドル（CUipcMemHandle）を開き、指定されたshape、strides、およびdtypeの配列として表すコンテキストマネージャー。stridesを省略したときは、1次元のC形連続配列であると見なされます。

戻り値デバイス配列を生成

コンテキストマネージャが終了すると、IPCハンドルは自動的に閉じられます。

CUDAアレイインターフェース（バージョン2）

cuda配列インターフェースは、さまざまなプロジェクトでのGPU配列のようなオブジェクトの異なる実装間の相互運用性のために作成されています。

このアイデアは、numpy配列インターフェイスから借用しています。

注意：現在、Python側のインターフェースのみを定義しています。

将来的には、コンパイルされたコードの情報を効率的に交換するために、

C側のインターフェースを追加する可能性があります。

Pythonインターフェース仕様

注意：仕様は変更になる場合があります。

__cuda_array_interface__アトリビュートは、次のエントリを含む必要がある辞書（dict型）を返します。

shape: (int, ...)

各次元のサイズを表すint（またはlong）のタプル。

typestr: str

NumPy配列インターフェースのtypestrと同じ定義です。

data: (int, boolean)

データは2タプルです。最初の要素は、Python int（またはlong）としてのデータポインターです。データはデバイスからアクセスできる必要があります。サイズがゼロの配列の場合、ここでは0を使用します。 2番目の要素は、Pythonboolとしての読み取り専用フラグです。

インターフェイスのユーザーは同じコンテキストにいる場合とそうでない場合があるため、最も一般的なケースは、CUDAドライバーAPI（または同等のCUDAランタイムAPI）でCU_POINTER_ATTRIBUTE_DEVICE_POINTERとともにcuPointerGetAttributeを使用して、現在で使用可能なデバイスポインターを取得することです。アクティブなコンテキスト。

version: int

エクスポートされるインターフェースのバージョンの整数。現在のバージョンは2です。

オプションのエントリは次のとおりです。

strides: None or (integer, ...)

ストライドが指定されていない場合、またはNoneの場合、配列はC形連続レイアウトになります。それ以外の場合は、各次元で次の要素にアクセスするためにスキップするバイト数を表す、int（またはlong）のタプルが明示的に指定されます。

descr

これは、より複雑な型を説明するためのものです。

numpy配列インターフェースと同じ仕様に従います。

mask

None、または__cuda_array_interface__を公開しているオブジェクト

Noneの場合、データ内のすべての値が有効です。マスク配列のすべての要素は、この配列のどの要素が有効であるかを示す trueまたは nottrue としてのみ解釈する必要があります。

これは、numpy配列インターフェースのマスクと同じ定義です。